## To create and send PBS files to run natural waterbody simulations ###

setwd("~/final_natural_water_files/PBS_Files")

script.file<-'~/evo-dispersal/art_wbdies/Pilbara/Pilbara_Natural.R'

for (ss in 1:100) {
  
    
    	fid<-ss
	    ##create the sh file
	    zz = file(paste("Pilb_Nat", fid,'.sh',sep=''),'w')
  	  cat('##################################\n',file=zz)
	    cat('#!/bin/sh\n',file=zz)
	    cat('cd $PBS_O_WORKDIR\n',file=zz)
  	  cat("R CMD BATCH --no-save --no-restore '--args File.ID=", fid, "' ", sep="", file=zz)
	    cat(script.file, " ", paste("Pilb_Nat", fid,'.Rout',sep=''), "\n", sep="",file=zz)
	    cat('##################################\n',file=zz)
  	  close(zz)
			
    	#submit the job
	    system(paste("qsub -m n Pilb_Nat", fid,".sh",sep=""))
  
}

########################
# And then concatenate the output
setwd("~/final_natural_water_files/Pilbara")
load(file="Natural waterbodies table.RData")
concat.times<-c()
concat.pops<-c()
for (ii in 1:100){
  load(file=paste("times colonised natural only", ii, ".RData", sep=""))
  concat.times<-c(concat.times, times)
  concat.pops<-rbind(concat.pops, points.colonised)
}

# proportion that colonised
sum(concat.times<100)/1000

#aggregate and calculate proportion of times that each point was colonised
nat.table<-cbind(nat.table, prop.col=tapply(concat.pops[,"Pres"], concat.pops[,"ID"], sum)/1000)
save(nat.table, file="Natural waterbodies table.RData")
colrs<-rev(heat.colors(10))
png(file="Nat waterbodies prop colonised.png", width=16, height=16, res=300, units="cm")
plot(nat.table[,"X"], nat.table[,"Y"], pch=19, bty="l", xlab="Easting", 
     ylab="Northing", col=colrs[cut(nat.table[,"prop.col"], breaks=10)])
legend("topleft", legend=c(0, 0.5, 1), fill=colrs[cut(c(0,0.5,1), breaks=10)])
dev.off()





